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We have developed a numerical code to study the evolution of self-gravitating matter in dynamic 
black hole axisymmetric spacetimes in general relativity. The matter fields are evolved with a high- 
resolution shock-capturing scheme that uses the characteristic information of the general relativistic 
hydrodynamic equations to build up a linearized Riemann solver. The spacetime is evolved with 
an axisymmetric ADM code designed to evolve a wormhole in full general relativity. We discuss 
• the numerical and algorithmic issues related to the effective coupling of the hydrodynamical and 

spacetime pieces of the code, as well as the numerical methods and gauge conditions we use to 
evolve such spacetimes. The code has been put through a series of tests that verify that it functions 
^ ' correctly. Particularly, we develop and describe a new set of testbed calculations and techniques 

designed to handle dynamically sliced, self-gravitating matter flows on black holes, and subject the 
code to these tests. We make some studies of the spherical and axisymmetric accretion onto a 
dynamic black hole, the fully dynamical evolution of imploding shells of dust with a black hole, the 
evolution of matter in rotating spacetimes, the gravitational radiation induced by the presence of 
the matter fields and the behavior of apparent horizons through the evolution. 



PACS numbers: 95.30.Sf, 95.30.Lz, 97.60.Lf, 04.25.Dm, 04.40.-b, 04.30.Db, 47.11.+j, 47.75.+f 
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OJQ. I. INTRODUCTION 

V~j ■ A. Overview 

An accurate description of the evolution of matter in a fully dynamical spacetime in complete generality is a 
longstanding and still unresolved problem in numerical relativistic astrophysics. With many important problems 
in urgent need of study, such as multidimensional core collapse, neutron star collisions or black hole formation, this 
outstanding problem represents a major gap in our ability to understand some of the most highly energetic processes in 
relativistic astrophysics. Numerous attempts have been made to develop such capabilities over the years, but only the 
spherical case can be considered essentially solved and even there it has not yet been widely applied. In higher 

dimensions, most of the studies have been restricted to the axisymmetric two dimensional (2D) case, and there much 
of the work has been devoted to hydrodynamical integrations on fixed general relativistic backgrounds |^-[f| ■ Even in 
axisymmetry, only a few attempts have been made to consider a fully self-gravitating, dynamic background [[l-I] |l7| . 
The three-dimensional case is also being studied, but the number of fully relativistic simulations is even smaller, 
with only a handful of fixed background hydrodynamical codes [ fbSlJlSfl and fully self-gravitating codes |2(J developed 
to date. In view of the difficulty of a fully relativistic treatment, approximations to a complete general relativistic 
approach are also being developed recently |^l|,^2| . 

The reasons for the reduced number of multidimensional (2D and 3D) computations are diverse, but we can 
enumerate the two most obvious ones: first, the inherent difficulties and complexities of the system of equations to 
integrate, the Einstein field equations coupled to the general relativistic hydrodynamic equations, which make this 
computation one of the most challenging ones in physics. Second, the immense computational resources needed to 
integrate the equations for 3D evolutions were not available in the past and are only starting to be achieved at present. 
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B. Axisymmetric General Relativistic Hydrodynamics 



In the axisymmetric case there exist a number of interesting astrophysical applications that can be addressed 
numerically, such as the rotational collapse of stellar cores in the supernova explosion scenario, the implosion of 
matter shells onto a black hole, the dynamics and stability of accretion disks or the fully dynamical spacetime version 
of the so-called relativistic Bondi-Hoyle accretion onto a moving black hole. Some of these problems have already 
tentatively been studied in the past within the framework of general relativity (see references below) although many 
scenarios still await for first (and detailed) computations. In this paper we take some early steps on that direction 
focusing on the methodology we will be using in the future for the study of those systems. 

In addition to the aforementioned scenarios, one of the most interesting systems to consider is the head-on collision of 
a neutron star with a black hole. The correct understanding of this simplified head-on model will undoubtedly serve as 
an important testbed for future three-dimensional codes. Hence, it is an important step that will aid in development 
of simulations of more complex scenarios, such as the coalescence and merging of spiral compact binaries. These 
catastrophic events are believed to be among the most promising sources of gravitational radiation to be detected 
by the gravitational wave interferometers to be operating around the turn of the century (LIGO, VIRGO, TAMA, 
GEO600). 

The use of general relativistic axisymmetric codes in numerical relativity has been largely devoted to the study of 
the gravitational collapse and bounce of a rotating stellar core and the subsequent emission of gravitational radiation. 
These investigations started with the work of Nakamura jL4|, who was the first to calculate a general relativistic 
rotating stellar collapse. In his calculation, he was able to track the evolution of matter and the formation of a 
black hole but the scheme was not accurate enough to compute the emitted gravitational radiation. Later, Stark 
and Piran 15 revisited this problem studying the collapse of rotating relativistic polytropic stars to black holes and 
succeeded in computing the gravitational radiation emission. Their code used the radial gauge and a mixture of 
polar and maximal slicing. The hydrodynamic equations were solved with standard finite difference methods with 
artificial viscosity p3 24 1. Evans [jl6| also studied the gravitational collapse problem but for non-rotating matter 
configurations. His numerical scheme to treat the matter fields was more sophisticated than previous approaches as it 
included monotonic upwind reconstruction procedures and flux limiters. Evans p6[ | was able to show that Newtonian 
gravity and the quadrupole formula for gravitational radiation were inadequate to study the problem. More recently, 
Abrahams et al. Jl7{ have solved the Einstein equations numerically for rotating spacetimes where the source of the 
gravitational field is a configuration of collisionless (dust) particles. They used the code to evaluate the stability 
of polytropic and toroidal star clusters. Because they did not have to take pressure forces into account, they could 
reduce the hydrodynamic computation to a straightforward integration of the geodesic equations. 

On the other hand, a few groups have concentrated on building codes to handle dynamical, axisymmetric gravi- 
tational fields in vacuum, without the complications of the matter fields. These have been applied notably to the 
collapse of gravitational waves to form a black hole [§5|,^6| and to the evolution of distorted, rotating, and colliding 
black holes fl27|-|30[| . But, even without the presence of matter fields, these calculations have proven very difficult for a 
number of reasons, not the least of which are the choice of appropriate gauge conditions and the presence of numerical 
instabilities encountered near the axis of symmetry. These are among the reasons that work in this area has been 
fairly limited. 



C. Our proposal 

We have embarked on a program to develop a series of fully relativistic codes to solve the coupled set of equations 
in multiple dimensions. This work builds on some years of experience in dynamic vacuum relativity (see, for exam- 
ple (3l|-|3^1) and in fixed background general relativistic hydrodynamics (see |]36| , ^2| , ^j37| ) . In spherical symmetry, 

] exploiting advances in hyperbolic treatments of both general relativistic hydrody- 
namics ]12J and the Einstein equations themselves |4(],|d|]. In the 3D case, a code is under development specifically 
for the study of neutron star collisions ]3l|f4^| , |43| ] . 

In the present paper we report first results obtained with a 2D axisymmetric numerical code designed to evolve 
rotating black hole spacetimes with self-gravitating matter. We believe that this code can have a number of interesting 
applications in the numerical modeling of accretion processes in black hole astrophysics. In addition, it will provide 
important testbeds for fully 3D codes now under development. This axisymmetric code is the result of the coupling 
of two previously developed production codes: an advanced general relativistic hydrodynamical code for stationary 
spacetimes and an axisymmetric black hole code to solve the dynamic Einstein field equations in vacuum. The 
integration of the different variables, spacetime and hydrodynamical quantities, is, hence, performed with a unique 
code. The two basic building blocks of the new code are extensively described in previous work |lj,^4|,|l],[52| . 



we have developed a code 
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The equations of general relativistic hydrodynamics are integrated with a modern high-resolution shock- capturing 
numerical scheme which relies on the knowledge of the characteristic information of the system of equations in order 
to build up a Riemann solver. These equations have been recently written as a hyperbolic system of conservation 
laws (with sources) for a general {3 + 1} metric Jl2| . The characteristic information of the system has been explicitly 
obtained. The Jacobian matrices are found to have real eigenvalues and a complete set of right-eigenvectors, therefore 
satisfying the definition of hyperbolici ty |4^| . The use of these methods, also referred to as Godunov-type methods 
(after the seminal idea by Godunov J46[) is becoming more important in recent years, due to a number of nice 
properties that other finite difference methods do not share (e.g., artificial viscosity methods). Those include a 
consistent treatment of discontinuous solutions (shock-capturing property) combined with being high-order methods 
in regions where the numerical solution is smooth. Recently, the use of these methods in relativistic hydrodynamics 
has permitted, for the first time, an accurate description of ultrarelativistic flows in different astrophysical scenarios 
(see, e.g., [ fL3"| , ]47[ | and references therein). 

The black hole piece of the code was originally developed by Brandt and Seidel |n],[52| and is based on the standard 
Arnowitt-Deser-Misner (ADM) [Q formulation of the Einstein equations as an initial value (or Cauchy) problem. 
The metric and extrinsic curvature components are evolved for the full set of the Einstein equations using a {3 + 1} 
explicit second order Runge-Kutta scheme with centered differencing. This code has a black hole built into the initial 
hypersurface of the spacetime. This avoids possible coordinate problems at the origin of the spherical coordinate 
system, since the black hole is constructed with an isometry that maps its interior, which contains a singularity at the 
origin, to the exterior, across a sphere at a finite radius. Hence, no reference need be made to the origin (as discussed 
in more detail below). This code can evolve a variety of spacetimes including rotating, distorted black holes. It also 
has a number of utilities built-in, such as routines to extract the various gravitational radiation modes and to track 
the motion of apparent and event horizons. 

The effective coupling of the two systems is through the source terms of the Einstein field equations. This allows us 
to integrate the whole system in a straightforward way - the metric and matter codes can simply take turns updating 
the variables. First, the hydrodynamical code takes a step, treating the spacetime metric as fixed and then the black 
hole code takes a step treating the matter fields as fixed. We regard this coupling of two mature codes as a starting 
point for algorithm development and forthcoming work in numerical relativistic astrophysics that will be explored in 
the near future. 

The organization of the paper is as follows: In section || we briefly review the ADM formulation of the space- 
time, discuss our particular choice of spacetime and hydrodynamical variables and describe the equations of general 
relativistic hydrodynamics. The initial value (Cauchy) problem with the presence of matter fields is discussed in sec- 
tion [II . We describe the details of the numerical code in section [V , where we address the specific issues concerning 



the coupling of the hydrodynamic equations and Einstein equations, as well as the gauge conditions used. Results 



and convergence tests of the simulations are presented in section VI . These include the spherical (Bondi) accretion 
of matter and the implosion (accretion) of a dust shell with the black hole. In the latter case, we also extract the 
waveforms of the gravitational waves induced by the presence of matter. Although in the present paper we mainly 
focus on the non-rotating case we will also present some results for rotating black hole spacetimes at the end of 
section Finally, section VII summarizes our main conclusions and outlines future directions of this work. 



II. PRELIMINARIES 
A. ADM formulation of spacetime 

We use the standard ADM (i^] formulation of the Einstein equations as the basis for our numerical code. Pertinent 
details of the formalism are summarized here, but we refer the reader to Ref. [ jf9| for a general treatment. A separate 
description of the theoretical details of both parts of the code can be found in Refs. |3lJ] and |l2|. Although for the 
sake of completeness we outline the most important points here, we refer the interested reader to these references as 
they will be extensively used throughout the present paper. 

In the ADM formulation, spacetime is foliated into a set of non-intersecting spacelike hypersurfaces. There are two 
kinematic variables which describe the evolution between these surfaces: the lapse function a, which describes the 
rate of advance of time along a timelike unit vector normal to a surface, and the spacelike shift vector /3 l that 
describes the motion of coordinates within a surface (throughout the paper Greek (Latin) indices run from to 3 (1 to 
3)). The choice of the lapse function and shift vector is essentially arbitrary (i.e., four degrees of freedom embodying 



the coordinate freedom of general relativity), and our choices will be described in section HE below. 
The line element is written as 

ds 2 = -(a 2 - Pi(?)dP + 2[3idx l dt + j ij dx i dx j , (I) 
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where 7^ is the 3-metric induced on each spacelike slice. Given a choice of lapse a and shift vector (3 i , the Einstein 
equations in the {3 + 1} formalism split into evolution equations for the 3-metric jij and constraint equations that 
must be satisfied on any time slice. The evolution equations are 

dtia = -laKit + Vifo + Vjfr (2) 
dtKti = -ViVj-a + a + K Ky - 2K lm Kf) 
+ f3 m V m K Z] + K m V 3 (3 m + K mj V % r 

- 8™ (Tij - ]pijT™ + ^PEjij^j (3) 

where Kij is the extrinsic curvature of the 3-dimensional time slice, Rij is the Ricci tensor of the induced 3-metric 
and Vj is the covariant three-space derivative. The matter terms involving Ty- and pe are defined below. 
The Hamiltonian constraint equation is 

R + K 2 - K l3 K l3 = 16np E (4) 

with R being the Ricci scalar. The three momentum constraint equations arc 

Vi (K ij - f ij K) = 8TrS j . (5) 

In the above equations, the quantities T 4 - 7 , S l and pE, the spatial components of the stress-energy tensor, the 
momenta and the total energy, respectively, are obtained by projecting the four stress-energy tensor using n^, the 
normal to the slice: 

-L* „ = g\ + n'n^ (6) 
T ij j_i v ± T n» (7) 

S i = - 1_ 4 *T» u n v (8) 
p E = 4 T^n^. (9) 

Note that Latin indices are raised and lowered with the induced 3-metric, e.g., K — K^^fi , = ftjij. As we 



discuss in section III, the constraint equations are used to obtain the initial data, and the evolution equations are 



used to advance the solution in time. 

B. Hydrodynamic Equations and Variables 

The matter fields appearing in the constraint and extrinsic curvature equations are computed via the local conser- 
vation laws of baryon number and energy-momentum, 

V p J 71 = (10) 

Vf/r^ = 0. (11) 

The rest-mass current and the perfect fluid stress-energy tensor of Eqs. ( |l0| ) and ( |ll| ) have the following definitions 

J" = pu" (12) 

and 

T" v = phu^v? +pg flu . (13) 

Here, p is the rest-mass density, p is the pressure, u M is the four-velocity of the fluid and h is the specific enthalpy, 
defined by h = 1 + e + p/p, where s is the specific internal energy. The spatial part of the fluid velocity is defined 
according to 



W^ + n") (14) 



such that v^n„ = 0. Given this condition, it follows that W is the Lorentz factor, W = —u^n,, = l/,/l — v„v^ and 
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± l X = Wv\ 
± l „ u^u" = W 2 v l v J . 



(15) 
(16) 



In Ref. ]l2| the general relativistic hydrodynamical equations are written, explicitly, as a hyperbolic system of balance 
laws in the framework of the ADM formulation. Starting from Eqs. ( |l0| ) and ( pi] ) and choosing an appropriate basis 

vectors for the spacetime , adapted to the Eulerian observers it is possible to cast the conservation equations into 
a more useful form fll2l: 



V, 



We can adopt a more compact notation as follows: 

V M f = s 

f = \j",Ti tu e^,Ti a 'n ll - J" 

Finally, this may be written in flux form by defining 

u = -f\/a 

giving us 



1 / d\/l u dV~9? 1 



-g \ dx° 



dx l 



— s. 



(17) 

(18) 
(19) 

(20) 
(21) 

(22) 



In the above equation 7 stands for the determinant of the 3-metric and yf—g = cx^Jj. Thus, the matter variables that 
are actually evolved in time are 



u= [D,S\r] 
= [J"n M) - V uT^n^T^n^ - J"n M ] 
= [pW, phW 2 v\ phW 2 -pW-P]. 



(23) 



For numerical applications we give the equations the form of a canonical balance law. Hence, we move quantities 
related to the spacetime metric to the RHS of Eq. ( p2| ) . Explicitly, we have 



du d(aP) 



Of 



dx l 



with 



s = as - 



u dy/j aP dy/^f 



(24) 



(25) 



C. Axisymmetric Coordinate System and Spacetime Variables 



Let us here introduce the notation for the spacetime variables, writing down the particular expressions used in 
axisymmetry from the general case considered in section IIA. Explicit details about this axisymmetric coordinate 
choice can be found in pl[ . The metric variables are given by: 



E sin 2 1 



:j 




c 


£sin 2 


B 


Fsin 


F sin 9 


Dsin 2 



lv<t> ^IH 

*' I C B Fswt-) J (2(m.) 
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and 




* 4 H c H B H F sm9 . (26b) 



In these expressions 77 is a logarithmic radial coordinate, and (8,4)) are the usual spherical angular coordinates. The 
relation between 77 and the standard radial coordinates used for Schwarzschild and Kerr black holes is given by 



r = e v (27) 



and 



/ M + a\ ( M-a\ 

where M is the mass of the black hole, a its angular momentum per unit mass, f a generalization of the Schwarzschild 
isotropic radius and r is the usual Boyer-Lindquist radial coordinate. This logarithmic coordinate allows one to 
impose outer boundary conditions at very large values of the standard radial coordinate. The conformal factor ^ is 



determined on the initial slice, as explained in section III below. Since we do not use it as a dynamical variable, it 
remains fixed in time afterwards. 



D. Boundary Conditions 

Conceptually, the computational domain consists of the region between two nested spheres, the throat of a black 
hole and a constant radius that is very far away (several hundred ADM masses). At the throat there is an isometry 
condition, which says that all variables (the metric and gauge variables, the matter fields, etc) at the interior of the 
throat can be calculated from their exterior values. 

Across the throat of the black hole, labeled by r\ = 0, we can demand the condition that the spacetime has the same 
geometry for positive 77 as for negative r\. This condition will build a wormhole into our spacetime. If we formulate 
the symmetry correctly, we will obtain simple boundary conditions for the throat that apply not only initially, but 
throughout the evolution. This boundary condition may be stated in one of two ways (in axisymmetry) to allow for 
different slicing conditions. Each choice must result in a symmetric and an anti-symmetric K n g to be consistent 
with the Kerr solution. Thus, the form of the isometry condition must be different for spacetimes sliced with a 
symmetric or an anti-symmetric lapse. For an anti-symmetric lapse the condition is 77 — > —77 while for a symmetric 
lapse the condition reads (77, 4>) — > (— 77, — 0). 

If a lapse that is anti-symmetric across the throat is desired, the metric elements with a single 77 index are anti- 
symmetric across the throat, while those with zero or two indices are symmetric. The extrinsic curvature components 
have the opposite symmetry of their corresponding metric elements. The shift /3 17 is anti-symmetric across the throat, 
while all other shifts are symmetric. 

The scalar matter fields, e.g., density or pressure, are simply symmetric across the throat, producing the same value 
at —77 as at +77. The momentum fields, S l are proportional to the shift divided by the lapse (/3 l jot) and therefore 
have the same symmetry as this quantity. 

If a symmetric lapse is desired, the metric elements with a single 77 index or single <j> index (but not both) will be 
anti-symmetric at the throat and all others will be symmetric. The extrinsic curvature components will have the same 
symmetries as their corresponding metric elements. The P v and (3^ shifts will be anti-symmetric, and the (3 9 shift 
will be symmetric. With these symmetries enforced, the initial data and all subsequent time slices will be isometric 
across the throat. One can verify that all Einstein equations respect these symmetries during the evolution if they 
are satisfied initially. 

In the ^-direction we can either evolve the whole region from 8 = to 8 = n or we can use an equatorial plane 
symmetry to increase the effective resolution of our simulation. 

All metric elements, extrinsic curvature components, shift components, and the S* momentum components, with a 
single 8 index are anti-symmetric across the symmetry axis. The remainder fields are symmetric. If we are not using 
an equatorial plane symmetry this condition applies both at the axis 8 — and at 8 = tt. 

At the equator (8 = 7r/2) there are two possible symmetries, the Kerr symmetry and the "cosmic screw" symmetry. 
For the Kerr symmetry, 8 —* tt — 8 , all metric components, extrinsic curvature components or shifts with a single 8 
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index are anti-symmetric. The remainder fields are symmetric. For the cosmic-screw type boundary conditions the 
symmetry at the equator is (#, 0) — > (7r — 9, — </>) and those metric elements, extrinsic curvature components, or shift 
components that contain one or one (f> index (but not both) are anti-symmetric. The remainder are symmetric. 
Finally, at the outer boundary a Robin condition is used for ^ H]. This condition gives the correct asymptotic 



behavior in the conformal factor to order r . For the metric given in the form (26a) ^ has the form 



* = e "/ 2 + (m/2)e- ri/2 + ... (29) 

and therefore obeys the differential equation ^ + 28^ = 2e* ? / 2 . The conformal factor is always symmetric at the 
throat, axis, and equator. 



E. Gauge Conditions 

Following Ref. J3l| , we will utilize the gauge freedom provided by the shift vector to reduce the number of spacetime 
variables that are evolved. Two of the shift components are used to eliminate the off-diagonal metric variable C and 
one shift component is used to eliminate the off-diagonal metric element-E. We note that this leaves one degree of 
gauge freedom unexploited. Also, due to the presence of both even- and odd-parity gravitational wave modes (or in 
other language, the "plus" and "cross" modes), the metric cannot be made completely diagonal. The 70^ component F 
carries information about the odd-parity wave modes that must be accounted for when rotational effects are included. 

As for the lapse choice, the code uses the time honored maximal slicing |5(J(that is, one of most commonly used 
slicing conditions in evolutions of black holes to date). The singularity avoiding properties of this slicing are charac- 
terized by the appearance of a limit surface at a distance from the black hole singularity that is dependent on the 



angular momentum |51 52], and, to a minor extent, on the form of gravity waves contained within the spacetime p4[ 
Maximal slicing is derived from the freely imposable condition that the trace of the extrinsic curvature should 
vanish throughout the evolution. We note that the Kerr solution in standard form is already maximally sliced with 
antisymmetric lapse, i.e. one which has the negative isometry sign across the isometry surface going into the black 
hole. Setting tnK = in the evolution for trif gives 

- V Q V a a + a (K ab K ab + 4tt (p E + S)) = 0. (30) 

This elliptic equation is solved numerically on each time slice during the evolution using a multigrid solver. This elliptic 
equation solver ]53[ ] is a semi-coarsening multigrid solver which does tridiagonal solves along the radial direction. It 
has proved quite robust and reliable in our numerical work to date. 

For the purposes of our numerical evolution, the previous considerations leave the metric variables A, B, D, F, 
and all six fTs as dynamical variables to be evolved. The various factors of s'md are included in the definitions to 
explicitly account for some of the behavior of the metric variables near the axis of symmetry and the equator. 

The condition on the shift used in our evolutions is, as in previous work on vacuum black hole spacetimes, 7 r) e = 
and 7^0 = 0. This choice simplifies the numerical equations and stabilizes the code. Since the Kerr shift allows the 
stationary rotating black hole metric to be manifestly time independent, one expects that for the dynamical case a 
similar shift will be helpful, and our procedure recovers the Kerr shift given the Kerr lapse, metric, and extrinsic 
curvature. 

We construct the shift condition by means of the evolution equations for the two metric variables we are setting to 
zero. Let us now consider how to implement the condition 7^0 = 0, 7,^ = - or, in terms of our variables, C = and 
E = 0. The relevant metric evolution equations are: 

d t C = 0= -2aH c + Ad n (3 9 + Bd B (i r > + Fsin<9<9„/3 (31a) 
d t E = 0= ~2aH E + Dd^ + Fd n f3 9 / sin 6. (31b) 

These equations can be combined to produce a single equation involving (3 V and f} e : 

2a^H c -^^-H E ^Ad,(3 e +^B-^pjd e f3^. (32) 

We can solve this equation by introducing an auxiliary function Q through the definitions: 

F> = d e n, (33a) 
f3 8 = d v n, (33b) 
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(following Ref. |f28|), producing an elliptic equation for the function fl: 

2a (h c - —j^—He\ = AtfSl + (b - ^\ dfa. (34) 

This equation is then solved by finite differencing using a numerical elliptic equation solver. The solution Q is then 
differentiated by centered derivatives to recover the shift components (3 V and (3 6 according to Eqs. (|33|). In practice, 
these shifts remain fairly small during the evolution. Their main function is to suppress the axis instability, as noted 
in Ref. |2S|] where a similar shift was used. 



Once ft is known, (3^ can be calculated by integrating Eq. (31b) 



^ = P(Kerr) + f d V (^H E - F8 v p e / sin 6) /D. (35) 

Only one boundary condition needs to be set (the outer boundary condition is most convenient) , and it is generally set 
equa l to t he Kerr value. The inner boundary condition, that fir must be symmetric across the throat, is guaranteed by 
Eq. (31b). This shift component is needed to keep the coordinates from becoming "tangled up" as they are dragged 
around by the rotating hole. Without such a shift the coordinates would rotate, leading to metric shear Q. This 
shift component, j3^, is typically larger than (3 V or (3 e . 

F. ADM Mass 

Within the axisymmetric coordinate system defined before, the ADM mass and angular momentum about the z-axis 
are defined to be |55[ 

M ADM =~i V a (Ve-^ 2 )dS a (36a) 
27r Js 

Pa = -^f s (H b a -l b a H)dS b . (36b) 

In terms of the variables defined in this paper these expressions yield 

Madm = - [ e ,l/2 (d v V - */2) sin ddO, (37a) 
Jo 

J = P<t> = \ f ^ H E\j^j- sin 3 Ode (37b) 



Because of this, the variable He is extremely important. It determines whether or not angular momentum is present 
in the spacetime. 

Although the ADM mass is defined strictly only at spatial infinity 1°, in practice we evaluate it at the edge of 
the spatial grid. As we use a logarithmic radial coordinate 77, this is in the asymptotic regime. While the angular 
momentum is, in principle, also measured at 7°, the presence of the azimuthal Killing vector makes it possible to 
evaluate J at any radius. However, unlike the vacuum case, the presence of matter does make it possible to transfer 
angular momentum through the motion of matter if there is angular momentum present in the matter initially. In all 
simulations in this paper, however, the angular momentum of the matter field is initially zero (i.e. S<f> — by choice; 
note that this does not mean that is zero, because that may include rotation of the coordinates). Therefore we 
expect that integration at any radius will yield J. We monitor this quantity during our simulations and use it as a 
test of the accuracy of our code. 



III. INITIAL VALUE PROBLEM WITH MATTER FIELDS 



Because the Einstein equations require initial values of the metric and extrinsic curvature which satisfy constraints 
(see Eqs. (||) and &), and because these constraints consist of four coupled non-linear partial differential equations, 
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obtaining a good starting point for an evolution requires a special technique - a conformal decomposition of the 
hydrodynamic and spacetime variables. 

The basics of this conformal decomposition of the initial data problem was first given by Lichnerowicz |5(| and 
later elaborated by Bowen and York |5j . Essentially, if one makes the following ansatz 

Ho = * 4 7ij (38a) 

PE = V- 8 pE (38b) 

Si = y- 6 S z (38c) 

Kij = ^- 2 H tJ (38d) 



and if a trace-free extrinsic curvature is used, then the Hamiltonian constraint decouples from the momentum con- 
straints, greatly simplifying the problem. 

Next, one normally assumes a form for the conformal metric, typically that it is flat, and then constructs a solution 
to the extrinsic curvature. Since this amounts to solving three equations for six unknowns it is customary to decompose 
the extrinsic curvature into a vector potential, Wj, (see |4t|) thereby reducing the number of unknowns to three: 

H i:i = 2V (i «7 i ) - 2/3V m w m 7!J . (39) 

Standard solutions to this equation with appropriate asymptotic behavior for linear or angular momentum have been 
obtained. See j57|] where solutions for single boosted or rotating black holes are described. 

Another simple solution to the momentum constraints can be found by simply inventing a form for w l , plugging it 
into the extrinsic curvature solution and then solving for the matter current term 



Si = ^ j H ij . (40) 



Because the momentum constraint equations are linear in Hij it is possible to combine them with the matter current 
solution obtained by inventing a w l simply by adding them together. One can add the form of the extrinsic curvature 
for a rotating black hole to the simple method described above, for example, to obtain a solution for initially flowing 
matter surrounding a rotating black hole. 

At this point we have nearly solved the problem. We have only two quantities left to determine: the conformal 
factor ^ and the conformal mass energy density pe- We shall consider pe first. Our primary consideration in choosing 
Pe will be to avoid unrealistic configurations of matter, i.e., matter which does not obey the energy condition 

S a S a < pi (41) 

and secondarily to find a solution which has matter flowing with a certain desired speed and along a chosen direction. 

For dust, the total velocity of the system can be calculated (thanks to the choice of the conformal factor made 
above) by the formula, 

?if = ?if = ViV i. (42) 

Pe Pe 

This is guaranteed to be less than one if we specify the conformal energy density in the following manner: 



PE = \JSi&+f E . (43) 

Here Je is an arbitrarily chosen function that respects the boundary conditions of our problem, and is everywhere 
greater than zero. The larger the value of /e the slower the matter will move. 

The situation is slightly different for the case of a perfect fluid. In this case we cannot know the precise value of the 
velocity until we have completely solved the problem. At this point the only unsolved quantity that remains is the 
conformal factor, 'J, and the only unsolved equation is the Hamiltonian constraint. We will solve the Hamiltonian 
constraint (a non-linear elliptic equation) for the conformal factor numerically using the multigrid solver mentioned 
above. Once this is done we can use Eq. (|SJ) and Eq. ( p"3| ) to obtain the following relations: 

S l S t = (phWf (W 2 - 1) (44) 
p E = phW 2 - P. (45) 
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We can then solve numerically these equations for p and W assuming a polytropic equation of state (i.e., constant 
spatial entropy condition) 

P = kp v (46) 
ph = p(l+e) + P (47) 
P = (r - l)pe, (48) 



with r being the (constant) adiabatic index of the fluid. It is assured by Eq. (44) that if a solution is found, ViV 1 < 1 
as in the case of dust considered above. 

IV. THE EVOLUTION CODE 

Now we turn to the description of the numerical code with which we solve the coupled set of equations presented 
in section |j. As already mentioned in the introduction, this code is the result of the merging of two previously 
existing independent codes. Each of these two matching pieces was originally developed to solve only one part of the 
problem. The final merged code is therefore capable of evolving either a vacuum spacetime, matter flows in a fixed 
background or a fully dynamic spacetime with evolving matter fields. Due to the ADM formulation of the Einstein 
equations, which we are adopting here for our numerical evolutions, the equations are a mixture of hyperbolic and 
elliptic equations. In consequence, it is not possible to write the full system as a single, unique system. Hence, the 
code evolves both fields in an alternate, almost independent, way. However, as we shall show below, we can still get 
second order convergence to the real solution using this approach. On the other hand, the alternate integration in 
time of both systems of equations allows one to use different numerical techniques for each one of them, choosing the 
best method for each piece. This is in fact what we do as we will explain later in this section. 

As stated previously, the hydrodynamical piece of the code makes use of a state-of-the-art modern high-resolution 
shock-capturing scheme. These methods are able to handle discontinuous solutions (e.g., shock waves) that could, 
eventually, develop in the flow, without using artificial viscosity to damp post-shock oscillations. To this aim th ey rely 



on the so-called approximate Riemann solvers. The code is capable of using either Roe's |5q| or Marquina's §59 60 
methods. In addition, the code can employ one of several different cell-reconstruction algorithms to increase the spatial 
accuracy of the hydrodynamic evolutions. Mathematically, these numerical schemes are built upon the characteristic 
information of the general relativistic hydrodynamics equations. Hence, the equations have to be written, explicitly, 
as a hyperbolic system of balance laws (as Eqs. ( p2| ) or (pi])). 

In the hydrodynamic integration of Eq. ( p4| ) the solution is updated in time, from time t n to time t n+1 , according 
to the following canonical conservative algorithm (written in ID for simplicity) 

= u " - + (49) 

where index j indicates cell centers and indices j ±1/2 indicate cell interfaces. In particular, and in order to increase 
the temporal order of the scheme, a high-order (typically second or third) monotonicity preserving Runge-Kutta 
method |nj is used to update the solution in time. For a third-order scheme the algorithm looks like this: 

u} = u" + AiL(u") (50a) 



3 3 ' ~"~V~.J 



u? = < + -AiL(uI?) + -AtZ,(uJ) (50b) 



" +1 = u? + ^AtL( Uj ") + iAtL(uj) + \ml{u)), (50c) 



with 



L = - (fj+1/2 A?' 1/2) +S ^ (51) 

A monotonic linear reconstruction of the cell centered values of the primitive variables provides second-order 
accuracy in space [ |62[ . Finally, the numerical fluxes across interfaces, fjii/a; are calculated, in the frame of the 
local characteristic approach, according to 
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fj±l/2 = "(^'±1/2 (wr) +f i± l/2(w L ) 



) (52) 



where wl and wr represent the values of the primitive variables at the left and right sides, respectively, of the 
corresponding interface. This state vector is defined as w = (p,Vi,s). In addition, {A n , r v n \n=x,. .5 are, respectively, 
the eigenvalues and right-eigenvectors of the Jacobian matrix of the system calculated at the interfaces from Wl 
and wr. Explicit general expressions can be found in ]l2|| . Finally, the quantities {Ajj n } n= i v .5, the jumps of the 
characteristic variables across each characteristic field, are obtained from 

5 

u(w R ) - u(w L ) = ^2 &v n r n . (53) 

n=l 

The tilde in some of the previous quantities indicates averaged quantities at the cell interfaces, algebraically computed 
from Wl and Wg. Further information can be obtained in p"2| ] and ply . 

For the spacetime part of the code we do not use a shock-capturing scheme. In our studies to date the metric 
variables are generally smooth so this presents no problem. However, in principle, any hydrodynamical shock can 
affect the spacetime part in terms containing first (time) derivatives of the extrinsic curvature components, as shown 
by Eq. (||). Hence, a shock-capturing scheme could still be a good choice. In addition, due to the slicing conditions 
and as an attempt to avoid any singularity appearing on the spacetime, the metric quantities can (and in fact do) 
develop large gradients (see, e.g., |3lf| ). In practice, these closely resemble the steep ones which characterize any 
real hydrodynamical shock. The use of shock-capturing schemes for the spacetime is, however, restricted to the 



recently developed hyperbolic formulations of the Einstein equations (see 64 6^,^] and references therein). We are 
also presently exploring the use of such advanced methods on the fully coupled equations in hyperbolic form |38) , |}9| , 
II- 

The spacetime variables are evolved for the full set of Einstein equations using an explicit second order Runge-Kutta 
scheme. Schematically the evolution looks like this: 

-n+i/2 = 7 « + 1a*7 (K n , 7 ") (54a) 

= + ^ ( R n ^ y, ) ( 54b ) 

= 7 " + A< 7 (iT l+1/2 , 7™ +1 / 2 ) (54c) 

K n+l =R n + M ft ^n+l/2 ) -n+lM ( 54d ) 



Spatial derivatives needed in the above equations are calculated using centered, second order finite differencing. This 
is a different scheme than we have employed in the past fl3l| . We switched the procedure in order to simplify the 
alignment of time levels. This different evolution procedure has proved to be as stable and reliable in our simulations 
as our previous scheme. 



V. ANALYSIS TOOLS 



Now we discuss the use of several tools we have developed to extract the physics from a fully relativistic simulation of 
dynamic spacetimes involving black holes. Unlike the standard cases involving hydrodynamics, where density, pressure 
and fluid velocities have well defined meanings, in a dynamic spacetime, the metric functions themselves, which are 
the functions actually evolved, do not carry direct physical or geometric meaning. Rather, they are related to the 
coordinate system and gauge in which the simulations are carried out. Therefore, physically relevant information 
about important quantities such as gravitational radiation or the mass of the black hole must be derived from the 
metric functions. 



A. Waveform extraction 



The gravitational radiation emitted is one of the most important quantities of interest in many astrophysical pro- 
cesses. However, the radiation is generated in regions of strong and dynamic gravitational fields, and then propagated 
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to regions far away where it can be detected. We take the approach of computing the evolution of the fields in a 
fully nonlinear way, while analyzing the radiation emitted in the regions where the system can be considered as a 
perturbation propagating on a fixed background. 

Under these conditions, one can appeal to the well developed theory of black hole perturbations. In this case one 
identifies certain perturbed metric quantities that evolve according to wave equations on the black hole background. 
However, the perturbed metric functions are also dependent of the gauge in which they are computed. Fortunately, 
there is a gauge-invariant prescription for isolating wave modes on black hole backgrounds, developed first by Mon- 
crief |]66| . The basic idea is that although the first order perturbed metric functions transform under first order 
coordinate transformations (gauge transformations), one can identify certain linear combinations of these functions 
that are invariant under such transformations. These gauge-invariant functions are clearly more directly related to 
true physics (which does not depend on coordinate systems), and in fact these functions obey the wave equations 
describing gravitational wave modes in linear theory. There are two independent wave modes, even- and odd-parity, 
corresponding to the two degrees of fr eedom, or polarization modes, of the waves. The metric used in this work allows 
both modes, as discussed in Sec. [IE above. 

In principle, one can consider a full set of gauge-invariants including both matter and gravitational wave fields, 
but our interest here is in evolving the hydrodynamics and spacetime within the fully nonlinear theory and simply 
extracting the gravitational wave information assuming a vacuum region away from the black hole. A waveform 
extraction procedure has been developed that allows one to process the metric and to identify the wave modes. The 
gravitational wave function (often called the "Zerilli function" for even-parity or the "Regge- Wheeler function" for 
odd-parity) can be computed by writing the metric as the sum of a background black hole part and a perturbation: 

o 

gap =9 a p ~\-h a p, where the perturbation h a p is expanded in spherical harmonics and their tensor generalizations. To 
compute the elements of h a p in a numerical simulation, one integrates the numerically evolved metric components 
g a p against appropriate spherical harmonics over a coordinate 2-sphere surrounding the black hole. The resulting 
functions can then be combined in a gauge- invariant way, following the prescription given by Moncrief f36f , leading 
to the gravitational wave functions denoted by ip. Then for each I — m mode, one can extract the waveforms of 
gravitational waves as they pass a "detector" at some fixed radius. This procedure has been described in detail for 
the case in this paper in [E3J31 3pJ, and more generally in Refs. p7|,|68[. 



B. Horizon tools 



We now turn to the topic of black hole horizons in numerical relativity and their application to the spacetimes 
considered in this paper. Horizons can be used in various ways to analyze the physics of a black hole system, and also 
as a check on the accuracy of the spacetime. We briefly discuss these issues below, and then apply horizon finders in 
testbed simulations below. 

There are two definitions of black hole horizons of interest to us: (a) The event horizon (EH) is the most commonly 
encountered term, defined loosely as the (closed, 2D) surface that separates those light rays that can escape the black 
hole's gravitational pull from those light rays that cannot. Exactly on this critical surface are light rays that never 
fall in to the black hole, and never escape to infinity. This surface is impossible to find on a given time slice, as 
photons that appear to be propagating (expanding) away from the black hole at one time may later find themselves 
falling back into the black hole if more mass-energy falls in, increasing its gravitational pull. Thus, the EH is generally 
an expanding surface composed of photons that will eventually find themselves trapped. Hence, locating the event 
horizon requires an entire evolution of a black hole spacetime. Methods do exist for finding event horizons [p9H72| and 
can be applied to spacetimes like those considered here, but we shall employ that application in future work, (b) The 
apparent horizon (AH) is defined loosely as the (outermost, closed, 2D) surface on which all outgoing photons normal 
to the surface are not instantaneously expanding away or converging towards each other: they have zero expansion. 
In stationary black hole systems, where no mass-energy will fall into the black hole, the AH and EH coincide, but 
generally the AH lies inside the EH. It is a convenient surface to locate on a given time slice, since one only needs to 
find a closed surface such that the expansion of all outgoing photon bundles have zero expansion. 

The expansion of a congruence of null rays moving in the outward normal direction to a closed surface can be 
shown to be [fT3| 

6 = Vi s l + Kij s l s j - tvK, (55) 

where Vi is the covariant derivative associated with the 3-metric 7^, s l is the normal vector to the surface, is the 
extrinsic curvature of the time slice, and trK is its trace. An AH is then the outermost surface such that 

9 = 0. (56) 
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This equation is not affected by the presence of matter, since it is purely geometric in nature, and so we can use 
standard horizon finders without modification for our current non-vacuum spacetimes. 

The key is to find a surface with normal vector s 1 satisfying this equation. There are many methods designed 
by now to locate apparent horizons in 2D, and we shall use the horizon finder developed for the black hole code 
used in this work, described in Ref. 44 3l]|^]. We refer the reader to those papers for details of the algorithm. In 



those references it was shown how one can use the dynamics and geometry of the AH surface to study the physics of 
dynamic, rotating vacuum black hole spacetimes. 

For the purposes of this paper, the most important physical quantity to be used is the area A of the apparent 
horizon. From this we compute an effective mass Mah of the black hole, defined in axisymmetry |74[ ] by 

o A 4ttJ 2 , , 

Ml„ = - + — . (57) 

We know that as the black hole approaches a stationary state at late times, the apparent and event horizons will 
coincide, and in that case the mass of the black hole is rigorously defined by the above formula. On e can use this 
black hole mass to study energy accounting: the total ADM mass of the spacetime, defined by Eq. ( |36aD above, should 
be equal to the final black hole mass plus any mass energy carried away by gravitational waves to infinity. This gives 
a powerful check on the overall global accuracy and consistency of the code, testing several crucial and independent 
aspects of the code and physics extraction models. These were very powerful tests in the vacuum case Q,|l],^2). We 
shall apply such tests below and show to what extent they are useful in these spacetimes with matter flows onto black 
holes. 



VI. TESTS AND CALCULATIONS 



In this section we present a summary of consistency checks and testbed computations that the code has successfully 
passed. We have first verified that the code reproduces previous results when either the matter fields are set to 
zero, or the evolution of the spacetime is turned off. We start by describing convergence tests of the code, for some 
particular initial matter configurations, to see the order of the method. This is a necessary requirement as no exact 
solutions to compare with exist. We then move on to consider simulations of the spherical accretion of dust and 
perfect fluid matter onto a dynamic black hole. In this application we show the differences between integrating the 
hydrodynamical equations fully coupled to the spacetime or in the simplified case of a sliced spacetime which does not 
react to the presence of matter. As expected, for sufficiently low energy density initial distributions, the spacetime 
evolution is totally unaffected by the matter content. Finally, we explore a few axisymmetric spacetimes and make 
some preliminary explorations of matter accreting onto both rotating and nonrotating black holes. 



A. Convergence tests 

We have performed a series of convergence tests on our code. We measure convergence along the line 6 = it /A for 
the constraints and other functions. Because we do not have data placed along this value of 9 we interpolate it from 
our existing data using a third order interpolation scheme. The convergence rate of the Hamiltonian and momentum 
constraints is based on data at two resolutions and the assumption that the true value is zero, and for any other 
quantity the convergence is calculated by comparing results obtained at three resolutions (keeping the ratio Ayy/A^ 
fixed) in a similar manner as reported in pjj j. In all cases, we assume that our functions are converging according to 
the formula 

X n = XTrue + k(nAT]) a (58) 

where k is a constant and n is either 1, 2, or 4 depending on whether we are using high, medium, or low resolution. 
The constraints are thus measured by substituting the numerically computed constraint for x n , zero for XTrue & n d 
the (radial) grid resolution for A77. The convergence of the constraints is given by 



log 



/log 2. (59) 



When we are not converging toward a known answer, for example when testing the convergence of the values of the 
field variables for the matter and the spacetime metric, we use the standard three point convergence formula 
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log 



x 2 - Xl 



/log 2. 



(60) 



In all cases that we show in this paper, we find that the convergence rate a is approximately 2 for most quantities, 
i.e., the code is globally second order accurate. It is worth mentioning that the two starting codes from which the 
present code has been built upon, that is the matter and the vacuum spacetime codes, were also, originally, second 
order accurate. The coupling of them both into a unique code has not diminished the final convergence order. 

For a sample case, analyzed in detail in Section VI C 2 below, we show convergence numbers in Table [j] and two 
plots of the constraint violation (See Figs. [[[ ||). The convergence is measured after an evolution time of 6Af. The 
resolutions used for this test are 100 x 10, 200 x 20, and 400 x 40 zones, in the radial and polar directions, respectively. 
The parameters describing the conformal density of the initial model arc (pb, k, p ma x, Vo, n ) — (0.01, 1, 3, 2, 2) and the 
angular momentum of the spacetime is 10 (to understand the meaning of these parameters we refer to Eq. (|6l|)). In 
the plots, we show the violation of the Hamiltonian constraint and the ^-component of the momentum constraint for 
the medium (xl/4) and high resolution runs. If the equations are converging at second order, then the two lines in 
each figure should be right on top of one another - as they are. 



B. Spherically Symmetric Simulations 

1. Dynamically sliced accretion 

We have first tested the hydrodynamical piece of the code against the analytic, spherically symmetric, Bondi 
accretion solutions ( |75| , see also Q for its relativistic extension) . We have verified that the code reproduces them, 
both for dust and a perfect fluid, when the background metric is kept fixed, to within a few percent for evolutions of 
1WM ADM . 

However, in dynamic spacetimes, one will not normally use a slicing that maintains a static metric, even if the 
underlying geometry has evolved to a stationary state. Therefore it is important to develop testbeds of known 
results, like the spherical Bondi accretion solution, but in non-analytic slicings which are commonly used in numerical 
relativity. We shall employ maximal slicing, but with different boundary conditions than those found in the stationary 
Schwarzschild metric. Therefore we compare the evolution of matter in a dynamically sliced spacetime which does not 
react to the presence of that matter distribution. In other words, the metric variables evolve completely independent 
of the hydrodynamical quantities (just as in vacuum spacetimes) whereas the latter evolve in a dynamic spacetime. 
This is a new type of testbed that we are developing for matter flows accreting onto dynamic, and dynamically sliced, 
black holes. 

We now wish to compare the analytic solution to the numerical one obtained in this dynamic slicing. It is nontrivial 
to do this in general, since one must compare invariant quantities at the same spacetime points in the two systems. 
Coordinate values will not suffice, since in the dynamically sliced case they are moving through spacetime. In the 
original slicing everything is a function of r only, while in the new slicing everything depends on space and time, 
i.e. on the new coordinates f and t. The rest-mass density, for example, is given by p{r) in the analytic solution, 
but in the new slicing it is an evolving quantity p(f , t) . Fortunately, it is possible to reconstruct r from f and t and 
thus to compare the time-evolved quantity, p{? , t) , with the analytic solution, p(r) . The value of r is simply the 
areal radius r — ^/jee- In Fig. || we see the result of an evolution with a non-Schwarzschild slicing. We plot the 
numerically evolved value of p (i.e. p(f, t)) with a dotted line. For reference, we also plot the initial value with a 
solid line. Note that when i = it is true that r — f. The line marked with circles (p eX act) represents the value of p 
as a function of the reconstructed r (i.e. p(r) — p(s/jee)). The solid line and the line for p eX act are really the same 
function, therefore, plotted with the areal radius as calculated at the times t = and t = 20M respectively. Clearly 
the numerically evolved solution (dotted line) agrees with the analytic (p eX act) solution to a high accuracy (and the 
constraints converged to second order). 



2. Dust accretion 



We now increase the complexity of the problem by "switching on" the matter fields appearing in the spacetime 
equations. We start showing the effects of a full coupling with the same matter fields of the previous section. Hence, 
initially, we assume a spherical distribution of dust with zero velocity and constant conformal energy density. Then, we 
solve the Hamiltonian and momentum constraint equations in order to obtain initial thermodynamic profiles satisfying 
Einstein equations. 
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First, we consider a small uniform density, p = 10 -2 , and solve the Hamiltonian constraint at t — 0, which, in turn, 
modifies the spacetime geometry. In Fig. [| we show the evolution of the radial metric function A up to a final time 
of 50Madm in intervals of 5M. We use a coarse radial grid of 150 zones and only one angular zone, with the outer 
boundary placed at n = 6.5. The solid line in this figure corresponds to the dynamically sliced Schwarzschild evolution 
discussed previously (a partly coupled evolution). Now we find, of course, a very different behavior from the standard 
Bondi accretion, where this metric function would remain time independent, which is merely a manifestation of the 
singularity avoiding maximal slicing used. The evolution of the radial metric function A for the fully coupled case 
is shown as a dotted line in Fig. ||. Note that for such a small density distribution, the effects of the matter on the 
spacetime geometry evolution are negligible. 

Now we show the effect of increasing the amount of dust. Fig. || corresponds to p = 1. Notice first that the peak 
of the metric function A at t = 50Madm (last curve) is located at a coordinate distance n slightly larger than 2 
regardless the value of p for the dynamically sliced Schwarzschild evolution (as one would expect). However, in the 
fully coupled evolutions, matter falls onto the black hole increasing its radius and causing the lapse to collapse more 
quickly in a larger interior region. As a result of this, the grid stretching effect which creates the peak in A occurs, 
in the high density run, at a more distant part of the grid than for the low density case. We have also carried out 
simulations with p = 10~ 3 . In this very low density case, as expected, the dotted and solid lines, or in other words, 
the fully and partly coupled evolutions, totally coincide. 

Let us now analyze the behavior of some characteristic relevant fields in these simulations. In the remaining of 
this section we consider the p = 1CP 2 case. Fig. || depicts the time evolution sequence of the lapse, a, the conformal 
radial metric variable A, and two hydrodynamical variables, the rest-mass density, p and the total velocity v. For 
this last quantity we plot its magnitude, multiplied by the sign of the radial component of the velocity (i.e. sign(v 7) ) 
y/v a v a ), in order to indicate if the matter is accreting onto the hole (negative values) or expanding away (positive 
values). The collapse of the lapse at the innermost zones freezes the evolution there at an early time, permitting us 
to avoid the singularity and continue the evolution until typical times of f=a IOOMadm- The radial metric function 
A shows its characteristic peak (a result of tidal forces on the grid points) which, at late times, leads to numerical 
problems in the integration of the metric part of the code. The density of the dust grows but rapidly settles in the 
interior region, due to the collapse of the lapse in that part of the grid and, although it is not visible in the plot, it 
continuously but slowly evolves in the outer zones. The evolution of the total velocity clearly shows that, initially, 
the matter accelerates rapidly, reaching ss 0.8c within the first 20M of evolution. However, due to the freezing of the 
lapse, this evolution is slowed and we never see it go past 0.9c even after 100A/ of evolution. If the evolution were 
continued long enough we would expect it to asymptotically approach c, which is the free-fall velocity at the horizon 
measured by an Eulerian observer in this particular coordinate system. 

In Fig. [?] we plot the time evolution of the apparent horizon mass for three different radial resolutions (150, 300 and 
600 zones). Correspondingly, in Fig. || we plot the location of the apparent horizon. We only show its location for the 
most resolved run (600 radial zones) as the results are almost identical for the lower resolutions. Note how sensitive 
the calculation of the horizon mass is to knowing its precise location. At lower resolutions, 150 and 300 zones, one 
sees a rapid, but spurious, growth of the horizon due to violation of the Hamiltonian constraint which actually causes 
it to go above the ADM mass within t = 50AI. In Fig. ^ we can see that initially 70% of the mass energy is contained 
within the apparent horizon, but by t = 50M almost all of this has fallen onto the black hole. 



3. Perfect fluid accretion 



Now let us increase again the complexity of the equations including the corresponding pressure terms which were 
absent in the dust case. Comparing to the accretion of pressureless matter, additional assumptions must be made in 
order to assign definite values to the density and pressure, as was described in section III. Furthermore, we make the 
choice that Pp^ 1 — constant initially. In Fig. [9] we plot the time evolution of the same quantities of Fig. ^| but now 
for the spherical accretion of a perfect fluid of adiabatic exponent 4/3. As before, we consider a constant conformal 
mass density p = 10~ 2 . Now, the role of the pressure makes the hydrodynamical evolution somewhat different from 
the pressureless case. This is clearly noticeable in the behavior of the velocity field. Initially all the material is at 
rest. The evolution of the velocity proceeds very rapidly and, in just 10M, we can identify part of the material falling 
towards the hole (negative values, using the same convention as for the dust case) and part of it expanding to infinity 
(positive values). The reason for this is the existence of a high pressure distribution of matter surrounding the hole 
and, hence, a pressure force that can, eventually, halt the gravitational collapse. The consequence is that part of the 
material bounces back towards larger values of n. This is in clear contrast with the dust case where no other force 
exists that can support gravity and the fluid can only freely-fall towards the hole. Finally, in Fig. [To] we plot the 
pressure and internal energy density evolution for this run. 
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C. Axisymmetric Simulations 

We present axisymmetric simulations showing the evolution of imploding (accreting) extended matter shells onto the 
black hole. We also discuss some results concerning the evolution of matter in rotating spacctimcs. The simulations 
we present here aim to show the performance and feasibility of our numerical tool. In future work we plan to carry 
over a detailed comparative and parametric study of the different astrophysical scenarios just outlined here. 



1. Imploding shells of dust 



As the first two-dimensional problem we study the evolution of an imploding shell of dust that radially falls to- 
wards the black hole. This problem has been studied semi-anal ytic ally in the past [fr7|, |Q where its astrophysical 
implications in gravitational wave astronomy were discussed. In [|77[ , the academic problem of computing the gravita- 
tional radiation from infinitesimally thin spheroidal sheets falling into black holes was discussed. In j7q], Shapiro and 
Wasserman extended these computations to finite size shells and pointed out the possible relevance of this process 
in non-spherical gravitational collapse to a black hole and in axisymmetric accretion onto a black hole. The most 
interesting conclusion of these semi-analytic studies was to demonstrate that the total energy emitted in gravitational 
radiation by a non-spherical dust cloud falling into a black hole is always less than that for a point particle of the 
same mass falling into the hole. This suppression was explained as a result from interference between waves emitted 
from different parts of the extended object. 

This problem has recently been studied numerically in ]79[| with a hybrid coupled code which employs linearized 



(perturbative) gravity (the Teukolsky equation with matter sources |80 ) and fully nonlinear hydrodynamics. The 
hydrodynamical piece of that code coincides with the one we use in the present fully non-linear code. In jT9| Pa- 
padopoulos and Font were able to demonstrate, in the linear regime, the suppression of gravitational radiation emission 
for a class of extended objects. Moreover, they showed that the radiated energy approaches an asymptotic value as 
the initial density distribution in the shell is made increasingly more compact. 

We present here the first numerical results of the fully non-linear evolution of an imploding shell of dust. Our aim 
is to capture the essential features of the problem. In order to do so we focus on a single initial model delaying for a 
forthcoming work a comparative study of the relevant parameter space of the problem (shell mass and compactness, 
shell-black hole separation, etc.). We are now mainly concerned in showing how the presence and evolution of the 
shell triggers the emission of gravitational radiation. 

We consider a Gaussian shell of dust whose conformal density distribution is parameterized by its location, 770, 
amplitude, p ma x and width, n, according to the following formula 



Pb 



~/W (e~^~^ 2 + e -«("+^) 2 ) sin" 9, (61) 



where pb is the background density. In particular we chose pb = 10~ 2 , p ma x = 10, « = 0.5, n = 4, and tjq — 2. We 
employ a grid of 300 x 25 zones in rj and 9, respectively. Then solve the Hamiltonian constraint, Eq. [I| The angular 
coordinate runs from to ir/2 and the radial one from to 7.5. The results of the simulation are plotted in Fig. |ll| 
for a final time of 100M. Here we plot a, A, p and v. One can again see the characteristic collapse of the lapse in 
the inner regions and the growth of the metric component A. The behavior of the total velocity is also similar to 
that found in the spherical accretion problem. The non-spherical aspect of this simulation is the initial distribution 
of the matter density. As this initial distribution has equatorial plane symmetry we only plot the results in the first 
quadrant. 



We plot in Fig. 12 profiles of the density along the axis (9 = 0) and for different times of the evolution. Clearly 
noticeable is the infall of the shell towards the hole accompanied by its collapse. As the shell does not have spherical 
symmetry, this implosion induces the emission of gravitational waves. We compute the even-parity £ — 2 and £ = 4 



modes of the emitted radiation and they are plotted in Figs. |13j and 14, respectively (odd-parity modes are absent in 
these non-rotating simulations). We also show in these plots a fit against the first two harmonics of the corresponding 
quasinormal modes. Clearly, they agree quite closely with the wave one expects to see for a black hole with a mass 
comparable to the entire dust-shell-plus-black-hole system. We emphasize that these quasinormal mode fits are not 
perturbative evolutions, as discussed above, but rather comparisons with known complex oscillation frequencies of 
black holes. The phase and amplitude of the two lowest modes are adjusted for a best fit to the obtained waveforms. 
We also compute the total energy radiated away by these two modes. For the £ — 2 we find 2.36 x 10~ 6 and for the 
I = 4 we get 6.35 x 10~ 9 . These values are normalized to the ADM mass of the spacetime. 

Next, we turn to another technique for studying a black hole accreting matter: examining the horizon dynamics. 
As discussed above, the black hole apparent horizon can give important information about the system. Its area is 
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related to the black hole mass, and combining the study of radiation emitted by gravitational waves can provide 
powerful checks on the overall energy accounting of the system. Energy conservation is very important in traditional 
hydrodynamics, but for dynamic spacetimes we must develop new techniques to account for radiation and the difficulty 
of localizing energy in general relativity. In this case one can use several completely independent measures of the 
energy in the system that should be related: the final mass Mah of the black hole, after all mass-energy has fallen 
in, will be given by the horizon area through Eq. (|57]) and the total energy radiated, E ra d can be computed from the 
Zerilli function ip. In principle these should add up to the total ADM mass of the spacetime, Madm, computed at 
large radius. This has been used extensively in vacuum. We now apply these techniques to black holes surrounded 
by matter. 

We consider a black hole plus matter spacetime with a non-flat background metric, hence inducing a more significant 
contribution of radiation energy. The specific non-flat background metric we employ is called the Brill wave metric, 
and it allows us to place a packet of gravitational wave energy of any desired width at any distance from the central 
black hole. For details of the construction of this metric see [Q. In Fig. [l5| we show various energy measures for 
this system. The thick lines shows the mass of the apparent horizon as a function of time. It grows as matter and 
gravitational wave energy fall in rapidly, and then settles down at roughly t = 20M. It is actually still growing slowly 
here because some matter is still falling in, but the amount is negligibly small. In order to compute the apparent 
horizon mass to the accuracy desired for this plot it was necessary to run at very high resolutions - in the final case 
shown the evolution was performed using 900 radial zones. As noted previously, this resolution is not needed to 
accurately compute the position of the apparent horizon, nor the energy in the radiation zone. 

We also show the ADM mass as a solid thin line, indicating the total mass of the spacetime. In principle, the black 
hole mass cannot exceed this limit, although it does at late times due to numerical error. Finally, we indicate the total 
radiated energy, computed through the Zerilli function, by a dotted thin line. The distance between the solid and 
dashed lines equals E ra d- If all mass energy has gone into the black hole, one should see Mah + E ra d — Madm- Fig. 
|l5| shows that this is quite closely achieved numerically, even though the total energy radiated is only 0.0874:M adm ■ 
The small gap between the energy radiated and the horizon mass is attributed to two effects: first, not all matter 
has actually fallen into the horizon by this time, and second, the apparent horizon mass will always be less than the 
event horizon mass in such cases. We regard this energy accounting issue as an important test and diagnostic of the 
physics of black hole accretion in dynamic black hole spacetimes. However, as one can see from Fig. |l5|, our studies 
indicate that this test is a very sensitive measure of global error. The crucial difficulty lies in resolving the peak in 
the metric function A that develops near the horizon. Small errors there translate into large deviations in the area 
calculation. Apparent horizon boundary conditions will aid this kind of study greatly. 



2. Rotating spacetimes 

Finally, we present some results concerning the evolution of matter in rotating spacetimes. The way in which these 
spacetimes are constructed and evolved is the same as in In a rotating spacetime, constant r\ observers along the 
equator are spinning around the black hole. In consequence, their fall through the horizon is slower and grid stretching 
is less near the equator. This is noticeable in the plots we show below. For the rotating cases below we generally 
choose a lapse which is symmetric across the throat, although in principle either choice (symmetric or antisymmetric) 
is possible. 

To show the behavior of the code when rotation is present we choose the same shell as in the previous section. The 
initial data set for the spacetime is based on a rotating Bowen and York |57|] black hole. We use this construction 
instead of Kerr for the sake of simplicity. We should note, however, that it contains radiation in the initial slice due to 
the construction of the initial data. However, the process of adding matter to a Kerr black hole and solving the initial 
constraints would also introduce radiation. Therefore, we choose the simpler initial data set provided by Bowen and 
York 0. 

Following jinj], we choose 

He = 3J, H F = 0, (62) 

where J = aM is the total angular momentum of the spacetime. We set initially J = 20. We use a grid of 300 
radial zones and 30 angular zones, with r\ max = 6.5. With the same shell parameters of the previous section we have 
an ADM mass of 6.87M which gives a maximum rotation parameter of a/M — 0.42. For comparison, we note that 
without the presence of matter this black hole system has an ADM mass of 4.75M. Hence, the surrounding matter has 
a significant effect on the black hole spacetime. However, the resulting black hole is not necessarily highly distorted, 
as it could result in a larger black hole with a small perturbation. 
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We evolve the initial data up to a time of 100M. By this time the final aspect of the lapse, metric function A, total 
velocity and density look very similar to the non-rotating run displayed in Fig. O. However, at much earlier times of 
the evolution, all variables present a more clear angular dependence as a consequence of the rotation. This is more 
pronounced as J increases. As an example we plot in Fig. |l6| the metric component A at t — 25Madm- 

The complexity of our system is now clearly greater than in previous models as we have additional quantities to 
evolve, e.g., the metric function F or the 4> component of the shift. The (f> component of the proper 4-velocity of the 
fluid is now different from zero as well. However, the corresponding 3-velocity component, as defined by 

air a 

should be zero, after being corrected by the shift. We find this to be the case in our numerical integration. 



We also use this run to test the ability of the code to conserve angular momentum as measured by Eq. 37b. We 
find that the initial angular momentum of the spacetime is accurately maintained to within 5% of the initial value, 
J = 20, at all points on the grid after t = 100M and converges away to second order. This behavior is indeed what 
one expects as there is no physical viscosity in the dust shell that could transfer the angular momentum. In fact, as 
stressed previously, the matter in this case is not carrying angular momentum, which can be seen from our initial 
data choice S$ — 0. The matter fields actually are measured by to the so-called Zero- Angular-Momentum-Observers 
(ZAMO's), which although rotating, do not carry angular momentum. 

We also compare the effects of the rotation and matter fields on the shape of the horizon. This is a further 
analysis technique that proved very useful in vacuum spacetime, which has also seen some use in matter spacetimes as 
well. Smarr ]81{] showed that for vacuum Kerr spacetimes, the horizon had an oblate shape parameterized by the spin 
parameter a. For larger a, the horizon becomes more oblate, as one might expect from naive considerations of spinning 
objects bulging at the equator. Previously, we found that for dynamic rotating black holes, the horizon oscillates 
about this oblate shape, settling down to its equilibrium value expected for a Kerr black hole of the appropriate 
mass and angular momentum. In fact, simpl y b y measuring the horizon shape, one could determine its mass, angular 
momentum, and oscillation frequency Q,|3l],|32| . We now apply this technique to black holes surrounded by accreting 
matter. 

For this purpose we plot in Fig. [l?] the ratio of the polar to the equatorial horizon circumference (C r — C p /C e ) for 
a sample of four runs: a vacuum run with J = 10, a run with low mass density and J = 10, a run with high mass 
density and J = 10, and a run with high mass density and J = 0. For each plot we also include a straight horizontal 
line corresponding to the value of C r for a Kerr black hole with the same value of cl/Madm- The low density matter 
distribution is given by Eq. ( |6l| ) with parameters pb = 10~ 2 , k = 1, rjo = 2 and p ma x = 1- Correspondingly, the high 
density matter distribution has the same values of K and 770 but the parameters pb and p m ax are ten times bigger. 

The point to notice in Fig. [l7] is that only the vacuum case settles down to the expected value of C' r for the 
Kerr spacetime given its angular momentum and ADM mass. The others settle down to something slightly different. 
In the low mass density case it is a little less spherical whereas in the high mass density case it is something a 
little more spherical. We expect that if all matter had fallen in the black hole, the standard Kerr result would be 
obtained. Clearly, the spacetime must be settling down to some quasi-stationary solution that corresponds to Kerr 
surrounded by matter. This is a very interesting point that should be explored further in future work. The effect of 
matter around a black hole on its geometry and oscillation structure has not received much attention, yet it could 
have important astrophysical consequences. As gravitational wave detectors begin to see waves from black holes, a 
particularly intriguing possibility is that they may carry information about not only the black holes themselves, but 
also about the astrophysical environment surrounding them |32| . 

The rotational implosion of the shell induces the emission of odd-parity gravitational waves in addition to even- 
parity modes. We plot in Figs. [l^ and ^0] the I — 2, 3 and 5 modes of the emitted radiation with a fit against the 
first two harmonics of the corresponding quasinormal modes. Again we note that these fits are made to the known 
quasinormal modes of vacuum black holes, although not all matter has crossed the horizon by this time. Further work 
should be done to study the effect of a "dirty" environment on the mode structure of black holes Q . 



VII. CONCLUSIONS 



We have presented a new numerical code to study the evolution of matter in black hole axisymmetric spacetimes in 
general relativity. Despite the well known "axis instability" of general relativistic axisymmetric codes we are able to 
evolve, during a reasonable amount of time into the future, different initial matter configurations. The two building 
blocks of the code, spacetime and hydrodynamics, are fully coupled through the source terms (and fluxes) of both 
systems. The extreme dynamical range of variation of the metric quantities, as shown typically in the peak of the radial 
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function A, the grid stretching or the formation of singularities, make the hydrodynamical integration quite a difficult 
enterprise. Within a dynamical spacetime framework, the computation becomes much more challenging than in the 
unrealistic case of an static background gravitational field, but it allows one to study the complete problem, including 
the effects of the matter on the black hole evolution and its corresponding structure and emission of gravitational 
waves. As this paper focuses on the development and testing of such a code and associated analysis tools to study 
the resulting physics, we defer detailed application to astrophysical scenarios to future investigations. 

The integration scheme used in the code is based on finite differencing the partial differential equations. For 
the hydrodynamic equations we have used an advanced high- resolution shock- capturing scheme built on approximate 
Riemann solvers. The integration of the ADM field equations was done with an standard explicit second order Runge- 
Kutta scheme with centered differencing. We have presented convergence tests of the code as well as a sufficient set 
of astrophysical applications. These include the spherical accretion of matter onto a black hole, the implosion of dust 
shells and the evolution of matter in a rotating black hole spacetime. We have also computed the waveforms induced 
by the presence of the matter in some of the aforementioned simulations. 

Because dynamic black holes accreting matter have not been studied previously, we developed a new series of 
testbeds appropriate for this problem and applied them to our code, (a) Building on the standard Bondi accretion 
on static black hole metric, which is an analytic solution, we showed how one can compare the numerical solution 
obtained by on a dynamically sliced background, (b) We also showed how one can computed radiation waveforms 
from the fully coupled matter-black hole system, which are emitted as the accretion induces oscillations in the black 
hole spacetime. (c) We applied a set of analysis tools developed to study the properties of black hole horizons in 
vacuum spacetimes to the accretion problem, and found them to be useful in studying the energy accounting of the 
entire black hole plus matter plus radiation system. We also studied the geometry of the black hole horizon as it is 
distorted by the presence of matter falling in. 

One particularly interesting point which emerges from these studies is the possibility that the matter surrounding 
the black hole perturbs it in a measurable way. The geometry of the hole is seen to be changed by the presence of 
matter, and it is possible that the radiation structure that one hopes to measure ultimately may be affected as well 
|S2f . This will need much more study in the future; a fully coupled hydrodynamics and spacetime code like the one 
developed here can address this problem in its full non-linearity. 

In subsequent work we plan to extend the results presented here performing detailed comparative and parametric 
studies of the different scenarios just outlined in the present investigation, including axisymmetric (non-spherical) 
accretion onto black holes and its effect on the structure of the black hole geometry and on the gravitational radiation 
emitted, the head-on collision of stars with black holes and a detailed comparison with recently developed perturbative 
treatments of matter flows around black holes |7jJ. We also plan to use these results and this code as a testbed for 
future computations with the three-dimensional coupled code, called Cactus, we are currently developing p3M3]. 
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var 


2-pt conv 


var 


3-pt conv 


Ham 


1.81 


P 


2.04 


Mom, 


1.97 


Srj 


1.76 


Mornj 


1.44 


A 


1.83 


Mom^ 


2.00 


B 


1.68 


trK 


2.04 


D 


1.89 



TABLE I. The result of a convergence test for the problem of an imploding dust shell with a rotating black hole with 
parameters (pi, k, pmax, Wo, n, J) = (0.01,1,3,2,2,10). The initial data is evolved for GMadm on grids of size 100 x 10, 
200 x 20, and 400 x 40 in 77 and 6, respectively. Convergence is nearly 2 for most quantities. The ^-momentum constraint is 
most affected by the axis instability and shows the lowest convergence. The number given is the average convergence value 
along the line 8 = n/4 for each quantity. 
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Hamiltonian Constraint 
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Position, T| 

FIG. 1. Numerical violation of the Hamiltonian constraint for the implosion of a dust shell onto a rotating black hole with 
J = 10, as measured along 6 — n/4. We plot the constraint violation at high resolution and 1/4 the constraint violation at 
medium resolution. For second order convergence these curves should lie on top of one another. 
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Momentum Constraint, r\ direction 



1 .0e-05 



-4.0e-05 



-9.0e-05 



-1.4e-04 




error, 400x40 
0.25*error, 200x20 



0.0 



2.0 



4.0 



6.0 



Position, r| 



FIG. 2. Numerical violation of the ^-component of the momentum constraint for the implosion of a dust shell onto a rotating 
black hole with J = 10, as measured along 6 = n/4. We plot the constraint violation at high resolution and 1/4 the constraint 
violation at medium resolution. For second order convergence these curves should lie on top of one another. 
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Bondi Solution with Alternate Slicing 

0.0040 p- 




0.0000 

10.0 20.0 

Schwarzschild 

FIG. 3. This figure shows the rest-mass density for a dynamically sliced spherical accretion of dust on a Schwarzschild 
background. The solid line represents the initial value - which is also the static analytic solution. The dotted line shows the 
density at a later time when the background spacetime is allowed to evolve with a dynamic slicing condition. The circles give 
the analytic density function but using ^fgee to calculate the Schwarzschild radius. The functions are identical to one part in 
1(T 5 . 
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FIG. 4. Time evolution of the metric component A for two different couplings between the hydrodynamics and metric 
components of the code. The solid lines indicate Schwarzschild (i.e., matter flowing on a dynamically sliced background which 
cannot be affected by the matter fields) and the dotted lines correspond to a ful coupling (i.e., fully self-gravitating matter 
that is allowed to affect the spacetime geometry). The last curve corresponds to t = 50Madm- The initial matter density is 
p — ICF 2 . Almost no differences are found between the two evolutions for such low density matter flows. 
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FIG. 5. Same as Fig. ^ but for an initial uniform density distribution p - 
different if the metric reacts or not to the presence of the matter. 
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1. Now, the evolution of the spacetime is noticeably 
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FIG. 6. Evolution of different metric and hydrodynamical quantities for the problem of spherical accretion of dust with 
p = 1CF 2 . The top panels show the lapse (left) and the metric component labelled A. The bottom plots show the rest-mass 
density (left) and the total velocity of the dust particles. 
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FIG. 7. Time evolution of the apparent horizon mass for the spherical accretion of dust problem. Each curve corresponds to 
a different radial resolution as indicated in the plot legend. 
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FIG. 8. Apparent horizon location as a function of time for the spherical dust accretion problem. 
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FIG. 9. Evolution of different metric and hydrodynamical quantities for the problem of spherical accretion of a perfect fluid. 
The top panels show the lapse (left) and the metric component labelled A. The bottom plots show the rest-mass density (left) 
and the total velocity of the fluid. 
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FIG. 11. Different metric and hydrodynamic quantities for the impact of an imploding shell of dust with the black hole. The 
solution is plotted at a final time of 100M. From top-left to bottom-right we show the density, total velocity, lapse and metric 
function A. 
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FIG. 12. Radial profiles of the density at different times of the evolution of an imploding shell. They correspond to an 
arbitrary constant value of the angular coordinate. Notice the collapse of the shell in progressive times. The thick solid line 
shows the initial profile. 
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1=2 Wave Mode Extraction 

Fit to Fundamental and First Harmonic of 1=7. 
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FIG. 13. This figure shows the numerically extracted i — 2 waveform (solid line) and the least squares fit to the two lowest 
i — 2 quasinormal mode for the imploding shell problem. 
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1=4 Wave Mode Extraction 
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FIG. 14. This figure shows the numerically extracted £ — 4 waveform (solid line) and the least squares fit to the two lowest 
£ — 4 quasinormal mode for the imploding shell problem. 
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Energy Accounting 
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FIG. 15. This figure shows the energy accounting for a simulation with a Brill wave and matter fields at a variety of 
resolutions. The Brill wave amplitude is 0.5, its width is 1.0 and its location is r/o = 2.5. The matter distribution is located at 
r)o = 2.0 and has an amplitude of 3.0. The background conformal density is 10~ 4 . 
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FIG. 16. Metric component A at t = 25Madm for a rotating black hole spacetime with J = 20. Notice the slightly less grid 
stretching near the equator as a consequence of the rotation. 
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FIG. 17. The ratio of the polar to equatorial circumference of the apparent horizon for a sample of rotating spacctimes. In 
each case, a horizontal line corresponding to the ratio for a Kerr black hole with the same ADM mass and angular momentum. 
The dotted line corresponds to a vacuum spacetime while the remaining three curves correspond to matter evolutions. One 
can easily see that only in the vacuum spacetime does the horizon settle down to something with the same shape as the Kerr 
black hole. 
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1=2 Wave Mode Extraction 



Fit to Fundamental and First Harmonic of 1=2. 
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FIG. 18. This figure shows the numerically extracted £ — 2 waveform (solid line) and the least squares fit to the two lowest 
£ — 2 quasinormal mode for the problem of the implosion of a dust shell onto a rotating black hole. 
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1=3 Wave Mode Extraction 



Fit to Fundamental and First Harmonic of 1=3 
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FIG. 19. This figure shows the numerically extracted i = 2 waveform (solid line) and the least squares fit to the two lowest 
= 3 quasinormal mode for the problem of the implosion of a dust shell onto a rotating black hole. 
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1=5 Wave Mode Extraction 



Fit to Fundamental and First Harmonic of 1=5 



0.00006 



0.00003 



B< 0.00000 



-0.00003 



-0.00006 



0.0 



Numerical Data 
Normal Mode Fit 




20.0 



40.0 60.0 
time (M) 



80.0 



100.0 



FIG. 20. This figure shows the numerically extracted £ — 2 waveform (solid line) and the least squares fit to the two lowest 
£ — 5 quasinormal mode for the problem of the implosion of a dust shell onto a rotating black hole. 



42 



